1 因果推断背景知识
我们知道，相关关系不等于因果关系，那么什么是因果关系，如何测度呢？
一般而言，因果关系是指变量 X 对变量 Y 在时间和空间上的实际影响。例如：广告促销活动对产品销售的影响有多大？
如果将广告促销作为一项干预（intervention），那么因果推断就是要准确评估该项干预对产品销量这一研究对象的因果影响。在因果推断研究中常遇到的一个问题是存在其他变量的混淆影响（confounding effect），以广告促销活动为例，可能是产品降价、宏观经济下滑或其他因素造成同一期间产品销量的变化，因此，怎么样将产品销量的变化归因为广告促销而非其他因素，就成为问题解决的关键。
通常，我们采取以下两种研究设计来估算某项干预对研究对象的真实因果影响：
1.自然/随机试验（Random-experiment）：这是推断干预的真实因果影响的最可靠方法，因为，用户随机直接触发的过程变化，可以进一步估算它对因变量的影响。但是，在很多情形下，这种随机试验或者无法实施，或者研究人员无法确保试验过程是完全随机发生。
2.准自然/随机试验（Quasi-experiment）：这种方法又称为计量因果推断法，主要通过应用统计计量手段控制混淆因素或应对反事实（counterfactual）难题，进而获得因果估算。这种方法又可以细分为常用的以下三种：
- 双重差分法（Difference in Differences，DID）或又称为”倍差法”
- 贝叶斯结构时间序列因果推断（Causal Impact Analysis based on BSTS ）
- 合成控制法（Synthetic Control Method, SCM）
这三种方法之间有何异同呢？理论和公式的讲解比较抽象、枯燥，我们通过一个分析实例来比较说明。
2 问题背景与数据集情况
2.1 问题背景
Abadie 和 Gardeazabal (2003) 最早提出了合成控制法，并发表论文以 1975年 西班牙 Basque Country（西班牙的一个自治区域） 发生的恐怖冲突事件对当地的经济影响评价为例来说明该方法如何使用。该文选择西班牙其他17个地区的经济数据作为恐怖事件发生前的合成控制变量，以使得合成的变量与 Basque Country 的经济变量相仿。
该数据集是面板数据，由18个分析单元构成，1个单元是干预单元（标号17，即Basque country），16个单元是控制（对照）单元（标号2-16,18），标号为1的单元是西班牙所有区域的平均值。数据字段包括：1个因变量，（人均 GDP：gdpcap）, 13个自变量（预测变量）,包括6个行业生产份额，5类识字或不同受教育程度人数，人口密度和投资比率等。地区编号和名称分别存放在字段 regionno、regionname 中，数据时间跨度42年（1955-1997）。
该数据集本来是用作合成控制法应用的示例，我们在此也将该数据集用于 DID 和 CI 方法的应用示例。
2.2 Synth包的安装
2.3 数据集情况
Basque 数据集情况如下：
- 数据记录的时间跨度 1955 - 1997
- 数据记录涉及 18 个西班牙地区和1个所有地区的平均值
- 恐怖事件发生的1975 年作为干预年
- 干预的地区是 Basque Country
- 经济影响的测度变量是人均 GDP （以千元计算）
数据集对应的数据字典（字段名称及含义）如下示意：
- regionno：Region Number
- egionname：Region Name
- year：Year
- gdpcap： real GDP per capita (in 1986 USD, thousands)
- sec.agriculture：production in agriculture, forestry, and fishing sector as a percentage of total production
- sec.energy：production in energy and water sector as a percentage of total production
- sec.industry：production in industrial sector as a percentage of total production
- sec.construction：production in construction and engineering sector as a percentage of total production
- sec.energy：production in marketable services sector as a percentage of total production
- sec.energy：production in Nonmarketable services sector as a percentage of total production
- school.illit：number of illiterate persons
- school.prim：number of persons with primary education or without studies
- school.med：number of persons with some high school education
- school.high：number of persons with high school degree
- school.post.high：number of persons with tertiary education
- popdens：population density (persons per square kilometer)
- invest: gross total investment as a share of GDP
# 清空环境
rm(list=ls())
begin<- Sys.time()
options(warn=-1) #消除警告信息
# 加载必要的包
library(Synth) # 做合成控制用
library(DT)
# 显示数据集
data(basque)
datatable(basque[,1:6], extensions = "Buttons",
options = list(dom = "Blfrtip",
buttons = c("copy","csv","excel")))3 三种估算方法的应用
3.1 双重差分法的应用
我们首先绘图看下人均 GDP 的变化情况，图 3.1：
# 加载必要的包
library(ggplot2) #加载 ggplot2 包
library(dplyr)
library(tidyr)
library(stargazer) # 输出统计表格
library(CausalImpact) # 做BSTS因果推断用
# 加载数据集
data(basque)
# 对数据集做清洗
unused <- c("sec.agriculture", "sec.energy" , "sec.industry" , "sec.construction" ,
"sec.services.venta" , "sec.services.nonventa", "school.illit", "school.prim",
"school.med", "school.high", "school.post.high", "popdens")
basq_clean <- basque[,!(names(basque) %in% unused)]
basq_clean <- basq_clean %>%
mutate(post = ifelse(year > 1975, 1, 0),
treat = ifelse(regionname == "Basque Country (Pais Vasco)", 1, 0),
regionname = as.factor(regionname)) %>%
dplyr::filter(regionno != 1)
# 输出图形
basq_fdid <- basq_clean %>%
dplyr::filter(treat == 1)
ggplot(basq_fdid, aes(x=year, y=gdpcap)) +
geom_line(color = "#F8766D") + theme_classic() +
geom_vline(xintercept=1975, color = "steelblue", linetype = "dashed") +
labs(title="",
y="人均GDP（千元）",x="时间（年份）", color = "Region") +
annotate("text", x = 1970, y = 9, label = "发生前", size =5, color = "steelblue") +
annotate("text", x = 1980, y = 9, label = "发生后", size =5, color = "steelblue") 
Figure 3.1: 1975 年恐怖事件发生前后 Basque 地区人均 GDP 变化趋势
从上图可以看到, 把恐怖事件作为一项干预，在干预时点前 Basque 地区的人均 GDP 基本是稳步上升，干预时点后，人均 GDP 先下降再回升。我们的研究就是要分析干预时点后的人均 GDP 下降幅度有多大是因为恐怖事件所造成的。
我们的分析思路如下：
首先，在原有数据字段基础上，增加2个字段：
字段1：post（区分干预时间点的虚拟变量：如果时间大于1975年，作为干预后阶段，标识为1，小于等于1975年，作为干预前阶段，标识为0）
字段2：treat（区分干预组和控制组的虚拟变量：如果该地区是 Basque Country ，那么就是干预组分析单元，标识为1，其他地区作为控制组分析单元，标识为0）
其次，用一阶差分来估算干预前后的人均 GDP 差异。为此，构建一阶差分方程，将人均 GDP 作为因变量，然后将虚拟变量 POST作为自变量，构建一元回归方程。
3.1.1 一阶差分估算
一阶差分估算的结果如下图示意：
# 计算一阶差分
f_did <- lm(data = basq_fdid, gdpcap ~ post)
stargazer(f_did, type="text")##
## ===============================================
## Dependent variable:
## ---------------------------
## gdpcap
## -----------------------------------------------
## post 2.484***
## (0.352)
##
## Constant 5.382***
## (0.252)
##
## -----------------------------------------------
## Observations 43
## R2 0.549
## Adjusted R2 0.538
## Residual Std. Error 1.153 (df = 41)
## F Statistic 49.921*** (df = 1; 41)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01从上表回归分析中可以看到：POST 虚拟变量的回归系数为2.484，说明干预后 Basque 地区人均 GDP 增长了2.484个单位。很明显，这不是我们需要的分析结果，因为我们需要的是干预发生后人均 GDP 下降的数值。
为什么会出现这种情况呢？仔细分析不难发现：Basque 地区人均 GDP 的变化受很多因素的影响，在我们重点关注的恐怖事件外，还可能存在如下一些混淆因素：
- 同期可能出台某个贸易法案，由此影响当地的商业和 GDP 值
- 当地军队的哗变
- 贪污腐败案爆发，政府部门停工等
解决办法是：将 Basque 地区的人均 GDP 变化与没有收到恐怖事件影响的其他地区做比对（这种研究设计考虑实际就是将受到干预影响的单元称为干预组，未受到干预影响的单元称为控制（对照）组），这一比对可以帮助我们剔除干预后的混淆因素，从而得到实际的因果影响，这就要应用双重差分方法来处理。
3.1.2 双重差分估算
双重差分方法应用的前提是：控制组的时间变化趋势与干预组未受到干预前的时间变化趋势相似。因此，干预前后，控制组与干预组对应曲线的斜率之差就是估算得到的干预影响（actual treatment effect）,这里的关键是假设控制组和干预组在干预之前的变化趋势一致。如果不符合这个前提，应用双重差分法就会存在问题。 在本例中，设置为控制组的区域是那些与 Basque country 比，历年的人均 GDP 方差最低的区域，这可以通过仔细查看控制组与干预组的人均 GDP 趋势曲线来获悉。本例中，Cataluna region 被选为最佳控制组所在地区，下图展示了控制组和干预组的人均 GDP 趋势。
#Picking the closest control group based on gdp
pre <- basq_clean %>%
dplyr::filter(post == 0) %>%
left_join(dplyr::select(basq_clean[basq_clean$post==0 & basq_clean$treat == 1, ], gdpcap, year),
by = c("year"= 'year')) %>%
mutate(perc_diff = (gdpcap.y - gdpcap.x) / gdpcap.y) %>%
group_by(regionname) %>%
summarise(gdp_var = abs(var(perc_diff))) %>%
arrange(gdp_var)
# Validating assumption
did_data <- basq_clean %>%
dplyr::filter(regionname %in% c("Basque Country (Pais Vasco)", "Cataluna"))
ggplot(did_data, aes(x=year, y=gdpcap, group = regionname)) +
geom_line(aes(color = regionname)) +
theme_classic() +
geom_vline(xintercept=1975, color = "steelblue", linetype = "dashed") +
labs(title="Basque 和 Cataluna 两地区 1955-1997 年人均 GDP 趋势",
y="人均 GDP（千元）",x="时间（年份）", color = "Region") +
scale_color_manual(labels = c("Basque (treated)", "Cataluna (control)"), values = c("#00BFC4", "#F8766D"))
上图可以看到：在恐怖事件发生前的1955-1975年期间，Cataluna 地区的人均 GDP 趋势与 Basque 地区的人均 GDP 趋势除了个别年份有所差异外，走势大体趋同。因此，将 Cataluna 地区作为控制组分析单元符合假设要求。
接下来，将人均 GDP 作为因变量，treat 和 post 的交互作为自变量，再做一次回归分析。这里的关键是应用 treat*post 交互关系来估算干预组干预后阶段的人均 GDP 与控制干预后阶段的人均 GDP 差值。回归分析的拟合结果如下表示意：
# Difference in differences
did <- lm(data = did_data, gdpcap ~ treat*post)
stargazer(did, type="text")##
## ===============================================
## Dependent variable:
## ---------------------------
## gdpcap
## -----------------------------------------------
## treat 0.139
## (0.376)
##
## post 3.339***
## (0.371)
##
## treat:post -0.855
## (0.525)
##
## Constant 5.244***
## (0.266)
##
## -----------------------------------------------
## Observations 86
## R2 0.607
## Adjusted R2 0.593
## Residual Std. Error 1.218 (df = 82)
## F Statistic 42.279*** (df = 3; 82)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01从上表的交互项系数（-0.855）看，在1975年的恐怖冲突事件发生后，Basque 地区的人均 GDP 下降了0.855个绝对值。不难发现，使用一阶差分估算和双重差分估算的结果差异明显。如果用一阶差分来估算干预影响将会造成严重错误。
双重差分方法有一定优势，也存在以下一些局限：
该方法基于传统的静态回归模型，假设因变量和自变量随时间不变，这明显不符合实际，因为大量的实际应用例子表明：回归模型中的系数应该是时变的（time-variant）；
大多数双重差分法将干预作为一个时点考虑，只分析干预前和干预后两个阶段的变化，但实际情况是，干预影响从产生到消失需要时间，干预影响本身不是一个时点，而是一个时间段。这样一来，双重差分法就无法应对了。
要解决上述问题，可以考虑使用时间序列中的时变因子模型，而贝叶斯结构时间序列（Bayesian Structural Time Series, BSTS）模型可以较好解决上述问题。下面我们通过应用了贝叶斯结构时间序列模型的 R 包— Causal Impact 来进行分析。
3.2 贝叶斯结构时间序列因果推断
贝叶斯结构时间序列因果推断是2015年由谷歌公司开发的一个前沿性因果推断方法，此处不再赘述方法的技术性原理，而是侧重讲解如何应用该方法。
采用贝叶斯结构时间序列进行因果推断的基本思路是：使用控制组的数据趋势来预测干预组在没有受到干预影响（这很明显是一个反事实）时的数据趋势，因此真实的因果估算就是干预组在干预影响后的实际值与干预影响后的反事实预测值之差。谷歌公司开发的因果推断 R 包— Causal Impact 使用贝叶斯结构时间序列模型来解释、预测干预组数据的反事实时空演化结果，本质上而言，Causal Impact 的思想非常近似于合成控制方法（Synthetic Control Method，SCM）。
在本例中，延续之前双重差分法的分析思路，控制组地区设定为 Cataluna 地区，将干预组和控制组地区对应的人均 GDP 值导入 Causal Impact R包即可输出得到结果。
# Causal Impact
basq_CI <- basq_clean %>%
dplyr::filter(regionname %in% c("Basque Country (Pais Vasco)", "Cataluna")) %>%
mutate(date = as.Date(paste0(year, "-01", "-01"), format = "%Y-%m-%d")) %>%
dplyr::select(date, regionname, gdpcap) %>%
spread(regionname, gdpcap)
names(basq_CI) <- c("date", "Basque", "another")
pre.period <- as.Date(c("1955-01-01", "1975-01-01"))
post.period <- as.Date(c("1976-01-01", "1997-01-01"))
impact <- CausalImpact(basq_CI, pre.period, post.period)
summary(impact) # 输出分析表## Posterior inference {CausalImpact}
##
## Average Cumulative
## Actual 7.9 173.1
## Prediction (s.d.) 8.6 (0.32) 190.3 (6.93)
## 95% CI [8, 9.3] [177, 203.9]
##
## Absolute effect (s.d.) -0.78 (0.32) -17.22 (6.93)
## 95% CI [-1.4, -0.18] [-30.9, -3.93]
##
## Relative effect (s.d.) -8.9% (3.3%) -8.9% (3.3%)
## 95% CI [-15%, -2.2%] [-15%, -2.2%]
##
## Posterior tail-area probability p: 0.00701
## Posterior prob. of a causal effect: 99.2986%
##
## For more details, type: summary(impact, "report")plot(impact) # 输出分析图
summary(impact,"report") # 输出分析报告## Analysis report {CausalImpact}
##
##
## During the post-intervention period, the response variable had an average value of approx. 7.87. By contrast, in the absence of an intervention, we would have expected an average response of 8.65. The 95% interval of this counterfactual prediction is [8.05, 9.27]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is -0.78 with a 95% interval of [-1.40, -0.18]. For a discussion of the significance of this effect, see below.
##
## Summing up the individual data points during the post-intervention period (which can only sometimes be meaningfully interpreted), the response variable had an overall value of 173.07. By contrast, had the intervention not taken place, we would have expected a sum of 190.28. The 95% interval of this prediction is [176.99, 203.92].
##
## The above results are given in terms of absolute numbers. In relative terms, the response variable showed a decrease of -9%. The 95% interval of this percentage is [-15%, -2%].
##
## This means that the negative effect observed during the intervention period is statistically significant. If the experimenter had expected a positive effect, it is recommended to double-check whether anomalies in the control variables may have caused an overly optimistic expectation of what should have happened in the response variable in the absence of the intervention.
##
## The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = 0.007). This means the causal effect can be considered statistically significant.上表中的绝对值即受干预（1975年恐怖事件）影响后的人均 GDP 与反事实人均 GDP 的差值（-0.76），意味着Basque 地区的人均 GDP 在恐怖事件发生后下降了0.76个单位的绝对值，相当于下降了 8.8% 。这与采用双重差分法计算的结果（-0.855）相近。
3.3 合成控制法的应用
3.3.1 方法简介
关于合成控制法的详细介绍见 Abadie and Gardeazabal (2003) 以及 Abadie, Diamond, and Hainmueller (2010, 2011, 2014). 该方法允许对单个分析单元（如企业、地区、国家等）受某个事件或干预影响进行评估分析，斯坦福大学著名学者 Athey and Imbens (2017) 称赞合成控制法为过去15年间最重要的政策评价模型。它提供了数据驱动的手段，通过设置干预组（受事件影响的分析单元）和控制/对照组（未受事件影响，但某些特征与干预组相似的分析单元），基于加权若干控制组与干预组的相似特性来建构合成控制单元，与干预组比较，若干控制组加权的结果要比单一控制组的结果更好。
与 Causal Impact 的估算类似，合成控制法也是要估算干预的实际影响。两个方法都要使用控制组来建构干预组的反事实（干预没有发生）结果。在本例中，干预组的反事实结果将通过控制组中的人均 GDP 以及其他相关协变量来完成。合成控制算法将控制组中的回归变量设置不同权重，以反映各变量对预测结果的不同影响，最终得到的实际因果影响就是干预后实际发生的人均 GDP 值与合成的反事实（干预没发生）人均 GDP值之差。
合成控制法与贝叶斯结构时间序列因果推断的区别在于：合成控制只使用干预前的预测变量来匹配，而 Causal Impact 使用干预前和干预后的预测变量时间序列来匹配，下图展示了数据集中 17 个地区人均 GDP 的变化趋势。
basq_synth <- basq_clean %>%
rename(Y = gdpcap) %>%
mutate(regionname = as.character(regionname))
ggplot(basq_synth, aes(x=year,y=Y,group=regionno)) +
geom_line(aes(color=as.factor(treat), size=as.factor(treat))) +
geom_vline(xintercept=1975,linetype="dashed", color = "steelblue") + theme_classic() +
labs(title="1975-1997 所有地区的人均 GDP 变化趋势",
y="人均 GDP （千元）",x="时间（年份）", color = "Treatment group") +
scale_color_manual(labels = c("Control", "Treated"), values = c("#F8766D", "#00BFC4")) +
scale_size_manual(values = c(0.5,1.5), guide = 'none')
图中可以看到,在干预前阶段，控制组所有地区的人均 GDP 走势和 Basque 地区类似，呈现上扬的态势。这说明基于控制组地区人均 GDP 走势来预测 Basque 地区人均 GDP 走势会相当准确。
3.3.2 应用步骤
R 中有好几个包可以实现从简单到很复杂的合成控制分析，以最为熟知的 Synth 包为例，合成控制法的应用有以下三步：
步骤1：应用合成控制法估计干预影响
首先明确干预前控制组参与合成控制的变量（这需要通过文献回顾和一点灵感找到），在本例中控制组每个地区有 13 个变量参与预测，如下：
1964-1969年地区总投资额/GDP （invest变量）;
1964-1969年工作适龄人口中不同教育程度（有教育经历:school.illit，有初等教育经历:school.prim，有部分高等教育经历:school.med，有完整高等教育经历:school.high，有研究生教育经历:school.post.high）的人口比例;
1961-1969年，6个工业部门占总生产产出的平均比例（以sec.前缀开始的变量）
1960-1969年，以1986年USD 不变价计算的人均 GPD 值（gdpcap变量）
1969年的地区人口密度，以每平方公里人数衡量（popdens变量）
合成控制法输出的结果如下图示意：
# synth
dataprep.out <-
dataprep(
foo = basque
,predictors= c("school.illit",
"school.prim",
"school.med",
"school.high",
"school.post.high"
,"invest"
)
,predictors.op = c("mean")
,dependent = c("gdpcap")
,unit.variable = c("regionno")
,time.variable = c("year")
,special.predictors = list(
list("gdpcap",1960:1969,c("mean")),
list("sec.agriculture",seq(1961,1969,2),c("mean")),
list("sec.energy",seq(1961,1969,2),c("mean")),
list("sec.industry",seq(1961,1969,2),c("mean")),
list("sec.construction",seq(1961,1969,2),c("mean")),
list("sec.services.venta",seq(1961,1969,2),c("mean")),
list("sec.services.nonventa",seq(1961,1969,2),c("mean")),
list("popdens",1969,c("mean")))
,treatment.identifier = 17
,controls.identifier = c(2:16,18)
,time.predictors.prior = c(1964:1969)
,time.optimize.ssr = c(1960:1969)
,unit.names.variable = c("regionname")
,time.plot = c(1955:1997)
)
synth.out = synth(dataprep.out)##
## X1, X0, Z1, Z0 all come directly from dataprep object.
##
##
## ****************
## searching for synthetic control unit
##
##
## ****************
## ****************
## ****************
##
## MSPE (LOSS V): 0.008864606
##
## solution.v:
## 0.02773094 1.194e-07 0.0000160609 0.0007163836 1.486e-07 0.002423908 0.0587055 0.2651997 0.02851006 0.291276 0.007994382 0.004053188 0.009398579 0.303975
##
## solution.w:
## 2.53e-08 4.63e-08 6.44e-08 2.81e-08 3.37e-08 4.844e-07 4.2e-08 4.69e-08 0.8508145 9.75e-08 3.2e-08 5.54e-08 0.1491843 4.86e-08 9.89e-08 1.162e-07# Two native plotting functions.
# Path.plot() plots the synthetic against the actual treated unit data.
path.plot(dataprep.res = dataprep.out, synth.res = synth.out,Xlab=" 时间（年份）",Ylab="人均 GDP（千元）")
abline(v=1975,lty=2,col="steelblue")
title("Basque 地区实际的人均 GDP 与合成的人均 GDP比较")
# Gaps.plot() shows the deviation between the synthetic and the actual over time.
gaps.plot(dataprep.res = dataprep.out, synth.res = synth.out,Xlab="时间（年份）",Ylab="人均 GDP 差异")
abline(v=1975,lty=2,col="red")
从图中可以看到，干预发生前，Basque 地区的人均 GDP 趋势曲线与合成的曲线吻合度很高，在干预发生后，人均 GDP 曲线逐渐偏离合成曲线，干预后的人均 GDP 曲线与合成曲线的差值就是干预的平均影响（Average Treated Effect, ATE）。干预后人均 GDP 曲线与合成曲线的均方根差异（Root Mean Squared Error，RMSE）值是 0.57 ,从合成控制法可以计算得到1975年恐怖冲突事件对 Basque 地区的影响是，人均 GDP 下降了0.57个单位值。
我们进一步考察17个地区中，哪些地区对合成控制变量的贡献权重大。
# 计算不同地区对合成变量贡献权重
synth.tables = synth.tab(dataprep.res = dataprep.out,
synth.res = synth.out)
# 有四个输出
# synth.tables$tab.pred 是比较干预前的预测变量值，合成控制变量值以及所有样本的均值
# synth.tables$tab.w 是干预前预测变量对合成变量的贡献
DT::datatable(synth.tables$tab.w)上表中可以看到，编号为10的 Cataluna 地区贡献了 85.1% , 编号为14的 Madrid 地区贡献了 14.9%，其他地区贡献为0，因此，使用合成控制法能通过算法来选择对合成变量贡献的协变量，不必像 DID 或 PSM 方法那样需要人为选定控制组单元。
步骤2：安慰剂检验（placebo test）
basque_dataprep<- function(treatment.id,controls.id){
dataprep.out <-
dataprep(foo = basque,
predictors = c("school.illit" , "school.prim" , "school.med" ,
"school.high" , "school.post.high" , "invest") ,
predictors.op = "mean" ,
time.predictors.prior = 1964:1969 ,
special.predictors = list(
list("gdpcap" , 1960:1969 , "mean"),
list("sec.agriculture" , seq(1961,1969,2), "mean"),
list("sec.energy" , seq(1961,1969,2), "mean"),
list("sec.industry" , seq(1961,1969,2), "mean"),
list("sec.construction" , seq(1961,1969,2), "mean"),
list("sec.services.venta" , seq(1961,1969,2), "mean"),
list("sec.services.nonventa" ,seq(1961,1969,2), "mean"),
list("popdens", 1969, "mean")
),
dependent = "gdpcap",
unit.variable = "regionno",
unit.names.variable = "regionname",
time.variable = "year",
treatment.identifier = treatment.id,
controls.identifier = controls.id,
time.optimize.ssr = 1960:1969,
time.plot = 1955:1997
)
dataprep.out$X1["school.high",] <-
dataprep.out$X1["school.high",] + dataprep.out$X1["school.post.high",]
dataprep.out$X1 <-
as.matrix(dataprep.out$X1[-which(rownames(dataprep.out$X1)=="school.post.high"),])
dataprep.out$X0["school.high",] <-
dataprep.out$X0["school.high",] + dataprep.out$X0["school.post.high",]
dataprep.out$X0 <-
dataprep.out$X0[-which(rownames(dataprep.out$X0)=="school.post.high"),]
lowest <- which(rownames(dataprep.out$X0)=="school.illit")
highest <- which(rownames(dataprep.out$X0)=="school.high")
dataprep.out$X1[lowest:highest,] <-
(100*dataprep.out$X1[lowest:highest,]) /
sum(dataprep.out$X1[lowest:highest,])
dataprep.out$X0[lowest:highest,] <-
100*scale(dataprep.out$X0[lowest:highest,],
center=FALSE,
scale=colSums(dataprep.out$X0[lowest:highest,])
)
return(dataprep.out)
}
dataprep.placebo <- basque_dataprep(10, c(2:9,11:16,18))
synth.placebo<- synth(
data.prep.obj = dataprep.placebo,
method = "BFGS"
)##
## X1, X0, Z1, Z0 all come directly from dataprep object.
##
##
## ****************
## searching for synthetic control unit
##
##
## ****************
## ****************
## ****************
##
## MSPE (LOSS V): 0.0003093978
##
## solution.v:
## 0.01405653 0.01342696 0.005037431 0.001792396 0.07939944 0.5166658 0.281674 0.08745683 9.4099e-06 0.0000164076 0.0000197882 0.0004392198 5.8367e-06
##
## solution.w:
## 0.0000109272 7.4858e-06 0.03591013 0.2715622 0.0000149049 0.2574583 3.0584e-06 2.667e-06 0.0000253151 2.2189e-06 4.9246e-06 0.4349737 0.0000181706 3.5542e-06 2.4455e-06# 以合成Basque中权重最高的地区作安慰剂检验
path.plot(dataprep.res = dataprep.placebo, synth.res = synth.placebo,Xlab=" 时间（年份）",Ylab="人均 GDP（千元）")
abline(v=1975,lty=2,col="steelblue")
title("合成Basque中权重最高的地区作为干预组，实际的人均 GDP 与合成的人均 GDP比较")
步骤3：排序检验（Permutation test）
# 将合成控制组各个地区作为干预组计算，检查合成控制组是否能够产生与干预组相似的结果
synth_results <- list()
dataprep_results <- list()
# 运行时间比较长
for (treatment_id in c(2:16,18)) {
output <- capture.output({
dataprep_result<-basque_dataprep(treatment_id, setdiff(2:18, treatment_id))
dataprep_results[[as.character(treatment_id)]] <- dataprep_result
synth_result<-synth(dataprep_result, method = "BFGS")
synth_results[[as.character(treatment_id)]] <- synth_result
})
}library(highcharter)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoochart <- highchart()
y_placebo<-matrix(NA,length(1955:1997),17)
# Basque地区的数据
x_Basque<-dataprep.out$tag$time.plot
y_Basque<-dataprep.out$Y1plot - (dataprep.out$Y0plot %*% synth.out$solution.w)
basque.mse <- apply(y_Basque^2,2,mean)
# 添加安慰剂试验数据
for (i in 1:length(dataprep_results)) {
if (i>=16)
treatment_id<-i+2
else
treatment_id<-i+1
x<-dataprep_results[[i]]$tag$time.plot
y<-dataprep_results[[i]]$Y1plot - (dataprep_results[[i]]$Y0plot %*% synth_results[[i]]$solution.w)
# 计算MSE
mse<-apply(y^2,2,mean)
print(mse)
# Exclude states with 5 times higher MSPE than basque
if (mse<0.5){
chart <- chart %>%
hc_add_series(data.frame(x = x, y = y), "line", hcaes(x = x, y = y), marker = list(enabled = FALSE),
name=paste0(treatment_id,"区域"), color="rgba(0,0,0,0.2)")
}
}## 2
## 0.1649906
## 3
## 0.03532337
## 4
## 0.3750022
## 5
## 1.693093
## 6
## 0.05041977
## 7
## 0.09071266
## 8
## 0.01189718
## 9
## 0.110995
## 10
## 0.2567429
## 11
## 0.009654997
## 12
## 0.5099131
## 13
## 0.005273292
## 14
## 0.6458686
## 15
## 0.07523803
## 16
## 0.07195226
## 18
## 0.02625146chart %>%
# 添加Basque地区的数据
hc_add_series(data.frame(x = x_Basque, y = y_Basque), "line", hcaes(x = x_Basque, y = y_Basque),
name="Basque地区", color="rgb(0,0,0)",marker = list(enabled = FALSE))%>%
# 标记处理时间为1975年
hc_xAxis(
plotLines = list(
list(
color = "red", # 虚线颜色为黑色
width = 2, # 虚线宽度为1像素
value = 1975, # 虚线在x轴上的值
dashStyle = "dash" # 虚线样式
)
)
)4 总结与启示
最后，我们将 DID、CI、SCM 三种方法计算的因果影响列示在下表中：
library(kableExtra)
labels <- c("双重差分法(DiD)", "贝叶斯结构时间序列因果推断(Causal Impact)", "合成控制法(SCM)")
values <- c(-0.85, -0.76, -0.57)
values_df <- data.frame(labels,values)
names(values_df) <- c("方法", "干预对人均 GDP 的影响估算")
kable(values_df)| 方法 | 干预对人均 GDP 的影响估算 |
|---|---|
| 双重差分法(DiD) | -0.85 |
| 贝叶斯结构时间序列因果推断(Causal Impact) | -0.76 |
| 合成控制法(SCM) | -0.57 |
可以看到，三种方法计算得到的干预影响值差别不大。当然，不存在给定唯一确定答案的干预影响及因果推断方法，实际应用时，我们要根据研究设计的限制，数据可获得性，以及干预影响的范围及复杂性等情况综合考虑。不过，我们可以综合应用这些方法，以确保研究结论的稳健性与可靠性。
本文的目的是讲解因果推断的三种方法(DID/CI/SCM)如何在 R 中应用，基于 Basque 数据集，为测算 1975年 Basque 地区发生的恐怖事件对该地区人均 GDP 的影响，我们分别采用了双重差分法（DID）、基于贝叶斯结构时间序列的因果推断法（Causal impact）以及合成控制法（SCM）来估算干预影响，通过上述分析，我们可以得到以下启示：
1.因果推断是一种十分重要且具有广泛用途的方法，它的一个目的是要定量评估某项干预事件/活动对研究对象的影响，因此，因果推断方法值得我们认真学习并熟练掌握；
2.因果推断方法是一种准自然实验方法，是对传统随机实验方法的重要补充，它的基本思路是将研究对象分为干预组和控制组，其中，干预组是受到干预影响的对象，控制组是未受到干预影响，但与干预组存在相似背景的对象，通过比较干预组和控制组在干预前后的变化，来估算干预事件/活动对干预的因果影响。
3.双重差分法、贝叶斯结构时间序列推断法、合成控制法有各自的优势和劣势，也有各自的适用条件，值得我们注意,列示如下：
双重差分法：假定控制组和干预组的数据时间趋势相同，将干预影响视作一个时点，而非一个时段。通过简单的回归，就可估算结果。如果数据样本量大，为确保干预组和控制组的数据样本特性相似，通常会先采用倾向得分匹配（Propensity score matching，PSM）法对样本进行匹配，然后再用双重差分来估算结果。有关倾向得分匹配法的细节与使用技巧可与老师课下交流，或看后续扩展版的讲解视频。
贝叶斯结构时间序列推断法：实质是合成控制法的一种，而且是较新的一种前沿方法。它的特点是考虑时变因素，即控制组和干预组数据随时间的变化，干预影响即可是一个时点，也可以是一个时段。而且，既可以将与干预组变量相关的其他变量纳入作为回归变量进行合成，也可以不考虑其他外生变量，只根据干预变量的历史事件序列来进行合成分析。这样就大大拓展了该方法在实际工作中的应用价值，因时间限制，该方法更多的细节与使用技巧可与老师课下交流，或看后续扩展版的讲解视频。
合成控制法：该方法近来在经济学实证、政治学实证研究中备受关注，常用于政策评估，事件评估等研究，因时间限制，本次讲解中未涉及安慰剂检验（ placebo test）的内容，更多的细节与使用技巧可与老师课下交流，或看后续扩展版的讲解视频。
本次数据分析耗时 1.34 分钟。